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Abstract 

We present a new scheme to extract numerically "optimal" interatomic po- 
tentials from large amounts of data produced by first-principles calculations. 
The method is based on fitting the potential to ab initio atomic forces of many 
atomic configurations, including surfaces, clusters, liquids and crystals at fi- 
nite temperature. The extensive data set overcomes the difficulties encoun- 
tered by traditional fitting approaches when using rich and complex analytic 
forms, allowing to construct potentials with a degree of accuracy compara- 
ble to that obtained by ab initio methods. A glue potential for aluminum 

obtained with this method is presented and discussed. 
PACS numbers: 34.20.Cf, 61.50.Lt, 64.70.Dv 
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While first-principles methods for computer simulation in condensed matter are rapidly 
improving in speed and accuracy, classical interatomic potentials continue to constitute the 
only way to perform molecular dynamics (MD) or Monte Carlo computations on systems 
with a very large size (number of atoms ~ 10^-10^) or for long simulation times {t ~ 
nanoseconds). With the advent of massively parallel machines and the proper computer 
codes, simulations on the mesoscopic scale appear feasible, allowing one to address a whole 
new range of problems in the physics of defects, surfaces, clusters, liquids and glasses. 
However, obtaining accurate and realistic potentials constitutes a challenging problem. It 
is now well recognized that fairly elaborate analytic expressions — involving for instance 
density-dependent terms, angular forces, or moment expansions — are necessary for a realistic 
description of most materials under different conditions (geometries, structures, thermody- 
namic phases). A typical potential is thus constituted by a number of functions combined in 
a complex way, and often nested one into another. Unfortunately, such powerful forms can 
make the task of fitting a potential to a given material quite formidable and cumbersome. 
There are often many possible ways to fit a set of experimental quantities within a given 
analytic framework, and rather arbitrary assumptions on the functions are usually made to 
reduce the number of parameters to a manageable level. Such assumptions could be the 
basic reason why potentials apparently good at T = sometimes fail at finite temperature, 
or for geometries or local conditions not considered when the fit was made. 

On the other hand, the development of first-principles methods — where forces on atoms 
are obtained by directly solving the electronic structure problem — has been very vigorous in 
the last decade , and evolved to a point where dynamical simulations of systems with 
of the order of 100-1000 and times t of the order of picoseconds are within reach for an ever 
increasing number of physical systems. Therefore, it seems compelling to construct a bridge 
between these two research lines, making use of the large amount of information that can 
be obtained by first-principles methods to construct reliable potentials for computations on 
a much larger scale. 

While one possible way to achieve this consists of trying to derive potentials from first- 
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principles theoretically by exploiting approximation schemes we proceed here along 

different, somewhat complementary lines. Namely, with realism of the final potential as the 
main goal, we present a new method to process a large amount of output of first-principles 
calculations (positions and forces), and combine this information with traditional fitting on 
experimental quantities, obtaining a potential by a numerical optimization procedure. 

We assume that the analytic form of the potential is defined by a number of single 
variable functions whose arguments are simply obtained from the atomic coordinates. Such 
is the case for essentially all the two-body and many-body potentials currently in use. For 
example, a glue potential 1^-^ 

^ ij i ^ j ^ 

is defined by a pair potential 0(r), an "atomic density" function p(r), and a "glue function" 
U{n). Other potentials may have functions of bond angles or of other quantities. Let 
a indicate the entire set of L parameters ai,...,aL used to characterize the functions. 
To determine the "optimal" set a* we try to match the forces supplied by first-principles 
calculations for a large set of different configurations with those predicted by the classical 
potential, by minimizing the objective function 

Z{a) = Zpia) + Zc{a) (2) 

with 

/ Af s -1 M Nk 

Z^(a)= S^iVfc EE|Ffc.(«)-FLI% (3) 

^ k=l ' k=l 1=1 
Nc 

Zc{a) = Y.Wr[Aria)-A:f. (4) 

r=l 

In Zp, M is the number of sets of atomic configurations available, is the number of atoms 
present in configuration k, Fki{a) is the force on the i-th atom in set k as obtained with 
parametrization a, and F^^ is the reference force from first-principles. Zc contains contri- 
butions from Nc additional constraints. Ar{a) are physical quantities as calculated with 
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parametrization a. A° are the corresponding reference quantities, which may be supphed 
either from the first principles calculation, or perhaps more likely directly from experimental 
data. Wr are weights which are chosen at convenience. 

The M configurations do not need to be related to each other, and in fact it is desirable to 
include input data relative to different geometries and physical situations (clusters, surfaces, 
bulk, solid, defects, liquid, etc.) in the attempt to achieve a good potential transferability. 
In practice, one can use samples from high temperature ah initio MD trajectories for various 
systems, thus obtaining a good representation of the regions of configuration space that they 
actually explore at finite T. 

Expression can be seen under two different points of view. If emphasis is given to Zp^ 
then the minimizing potential appears as the "best" approximation of the first-principles 
system, with Zc acting as a guide towards the correct region in the space of parameters. 
In fact, many properties (for instance, the formation energy of a defect) are not easily 
determined from the forces alone. Moreover, one could include terms aimed at correcting 
shortcomings of the first-principles method in use. If emphasis is given, instead, to Zc-, then 
the method looks like a conventional fit, but where the Zp term relieves the researcher from 
the enormous burden of guessing the shape of the functions constituting the model. ^4° and 
Wr can be maneuvered for tuning. 

It should also be noted that invariance properties of the Hamiltonian must be recognized 
and taken care of by imposing additional, dummy constraints. For instance, a glue potential 
(H) is invariant under the transformations (a) p(r) Ap{r),U{n) U{n/A), and (b) 
0(r) — > 0(r) + 2Bp{r),U{n) — > U{n) — Bn 0. The two constants A and B are arbitrary 
and must be fixed by additional conditions, that can be enforced as further quadratic terms 
in eq. (0), as if they were constraints for physical properties. In contrast with the latter, 
these terms exactly vanish at the minimum. 

In the present realization, the single variable functions constituting the potential are 
defined as third-order polynomials (cubic splines) connecting a set of points preserving 
continuity of the functions and of their first two derivatives across the junctions. The pa- 
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rameters ai are a one-to-one mapping to the points chosen on the basis of computational 
convenience P,p!0|. In the simplest case, (3i = ai. Particular boundary conditions, such as 



requiring a function and its first derivative to be zero at a cutoff distance Re-, are directly 
incorporated into the parametrization. A number of parameters of the order of 10-20 per 
function seems to give sufficient flexibility, and in fact using finer grids may give rise to 
noise problems (oscillations of the functions depending on the input set) in the spline-based 
formulation. 

The computational engine of the method is a multidimensional minimization procedure 
for the objective function @). To be prepared to deal with the presence of multiple local min- 
ima, we have implemented a simulated annealing algorithm in parameter space (described 
in [p|,p!0[). However, we found it to be necessary only when starting from an initial guess 
very far from the optimal one. The basin of attraction of the optimal potential is sufficiently 
broad that a simple quasi-Newton method is adequate to reminimize Z after small adjust- 
ments to the values of A° and Wr-, or changes in the set of first-principles configurations. In 
a typical minimization run, Z is evaluated a few thousand times. The force computations 
are carried out by using standard MD techniques such as neighbor lists to decrease computer 
time, and the computational resources can be compared to those requested by a classical 
MD code. In preliminary tests using MD trajectories generated by classical potentials [§, 
this scheme has proven to be able to reconstruct exactly the original potentials without any 
further assumption beyond the analytic form — within the precision allowed by the spacing 
between spline knots, and within the range of the function arguments sampled by the input 
data. 

The first application of the force-matching method presented here is for aluminum. The 
reasons for this choice are threefold: (i) due to the absence of d electrons, AI can be studied 
quite easily and accurately with present first-principles methods; (ii) the metallic character 
suggests the use of the relatively simple glue model (p, even if it has well known limitations 
(lack of angular forces); (in) in spite of the simpler electronic structure, glue-like potentials 
obtained so far for AI seem to be less accurate than those for noble metals, and the validity 



of glue schemes for Al has been recently questioned [TT|] , making it worthwhile to investigate 
this issue further. 

The potential has been parametrized by a total of L = 40 parameters (spline knots), of 
which 14 for (f){r) and p(r), and 12 for U{n). The functions reach zero at r = Rc, and at 
n = 0, by means of fixed additional spline knots. The first-principles data used as input have 
been extracted by trajectories of MD simulations using the local orbital density functional 
scheme described in 0. We have processed a total of M = 85 sets of atomic configurations, 
of which 7 represent a bulk system with a vacancy (A^ = 107) at T = 100 K, 10 the same 
system at 1750 K (undergoing melting), 20 an equilibrated bulk liquid (A^ = 108) at 2650 
K, 13 a (100) slab (A^ = 108, 100 K), 10 a A^ = 150 cluster at 1000 K, 25 the same cluster 
in the liquid state at 2200 K, for a total of 10633 force vectors included in Eq. (|]). 

We have used = 32 additional constraints, 8 of which for the cohesive energy, the 
equilibrium lattice spacing Oo, the (unrelaxed) vacancy formation energy, the (unrelaxed) 
(111) intrinsic stacking fault energy, the (unrelaxed) (111) surface energy, the bulk modulus 
and the shear moduli Cu — C12 and C44, 22 to fit the energy and pressure to the universal 
equation of state [O at 11 different lattice spacings (a/oo = 0.90, 0.94, 0.97, 1.05, 1.11, 1.20, 



1.30, 1.40, 1.50, 1.60, 1.75), and the remaining 2 are related to the invariance properties of the 
potential described above. The weights Wr assigned to the constraints and the cutoff radius 
Rc for (j){r) and p(r) have been adjusted by a trial-and-error procedure, where potentials 
were generated by minimizing (^) and then run through a test suite including evaluation of 
relaxed energies of defects and surfaces, surface relaxations, thermal expansion and a melting 
point estimate by zero pressure MD simulations. The final potential, shown in Fig. |l], has 
Rc = 5.56 A (between the 3rd and the 4th neighbor shell in the fee crystal), and corresponds 
to Z = Zp + Zc = 0.029 + 0.003 = 0.032 (eV/A)^. ^/Z^ - 0.17 eV/A is the root mean 
square (rms) deviation of force components, to be compared with the rms force component 
in the input data, 0.92 eV/A. Such error is at least one order of magnitude smaller than 
that typical of empirical models |TB|. Some properties of the potential are listed in Table |. 



It should be noted that no constraint was imposed on the phonon frequencies at the zone 



boundary, so that they are mostly determined by the force-matching term. Surface energies 
are somewhat lower than in experiment, although they are higher than those predicted by 
other potentials of the same class |p . Surface relaxations are very realistic — in particular, the 
rather uncommon surface expansion of Al(lll) is obtained — except that the contraction 
of Al(llO) is not as large as in experiment. This discrepancy is already present with the first- 



principle method in use and has been simply transmitted to the potential. The thermal 
expansion behavior, obtained by MD at zero pressure, is shown in Fig. ^. The melting 
temperature, corresponding to the discontinuity, has been determined with a precision of 
about 3 K by achieving solid-liquid coexistence in a constant enthalpy run for a system with 



10752 particles, following the technique described in [jT5[. Thermal and melting properties 
are in remarkable agreement with experiment. 

In conclusion, this first study shows that the force-matching method is a very effective 
tool to obtain realistic classical potentials with a high degree of transferability for systems 
which the ab initio calculation technology is capable of treating. The number of such sys- 
tems is rapidly increasing as electronic structure methods are improved and the computing 
power increases. The numerical optimization procedure at the heart of the method is ex- 
pected to be well suited to handle easily the rich and complex analytic forms — including 
angular-dependent terms — required for a realistic modelling of covalent bonds, and consid- 
ered difficult to fit so far. In fact, it could also be used to compare quantitatively different 
functional forms, on the basis of their accuracy in reproducing the ab initio forces. It 
would also be of considerable help in the fitting of alloys, where the number of experimental 
properties available is usually rather limited. 

We are indebted to Dave Drabold, Cathy Rohrer, Wei Xu and Sang Yang for useful 
discussions and suggestions. We especially thank Dave Drabold for providing the ab initio 
data. This work has been carried out under the U.S. Department of Energy Grant No. 
DOE-BES (0) 76ER01198. 
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FIGURES 

FIG. 1. The three functions constituting the optimized glue potential for Al. 

FIG. 2. Lattice parameter a as a function of temperature for our optimized potential (solid 
line), compared with experimental data (dotted line). The jumps indicate the volume change on 
melting. In the liquid region, a is defined as (40)^/"^ where Vl is the atomic volume, as in the fee 
crystal. 
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TABLES 

TABLE I. Experimental and calculated (with the optimized potential) values for equilibrium 

lattice spacing, cohesive energy, bulk modulus, elastic constants, phonon frequencies at the points 
X, L and K of the Brillouin zone, vacancy formation and migration energies, intrinsic (111) stacking 
fault energy, surface energy and surface relaxation between the two outmost layers for the (111), 
(100) and (110) surfaces, melting temperature, latent heat and volume change on melting. All the 
energies are at T = and include relaxation effects. 
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di2 (111) (%) 
di2 (100) (%) 
di2 (110) (%) 
(K) 

Lm (eV/atom) 

Ay„ (%) 



+0.9 ±0.7 § 
-1.2 ± 1.2 
-8.5 ± 1.0 ' 
933.6 
0.1085 
6.5 



+0.9 
-1.5 
-4.6 
939 + 3 
0.1054 + 0.0003 
8.4 + 0.1 



^Extrapolated classically to T = from data in Ref. [16 
'^Frequencies at 80 K from Ref. 



'^Ref. [118|. 
■^Ref. [|19|. 



^Ref. |20|. 

^Estimates for an "average" orientation, Ref. 



SRef. [|2|. 
^Ref. |23|. 

^Ref. H. Ref. l2|] reports -8.4 + 0.8. 
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